Systematic simulations of modified gravity: chameleon models 
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In this work we systematically study the linear and nonlinear structure formation in chameleon theories of 
modified gravity, using a generic parameterisation which describes a large class of models using only 4 param- 
eters. For this we have modified the A^-body simulation code ECOSMOG to perform a total of 65 simulations for 
different models and parameter values, including the default ACDM. These simulations enable us to explore a 
significant portion of the parameter space. We have studied the effects of modified gravity on the matter power 
spectrum and mass function, and found a rich and interesting phenomenology where the difference with the 
ACDM paradigm cannot be reproduced by a linear analysis even on scales as large as k ~ 0.05 /iMpc^^, since 
the latter incorrectly assumes that the modification of gravity depends only on the background matter density. 
Our results show that the chameleon screening mechanism is significantly more efficient than other mechanisms 
such as the dilaton and symmetron, especially in high-density regions and at early times, and can serve as a 
guidance to determine the parts of the chameleon parameter space which are cosmologically interesting and 
thus merit further studies in the future. 



I. INTRODUCTION 

Two plausible alternative explanations to the observed ac- 
celerating expansion of our Universe are dynamical dark en- 
ergy d] and modified gravity In both classes of theories, 
a scalar field has been used as the most common dynami- 
cal origin of the acceleration of the Universe. This, however, 
comes at a price: in many theories, especially modified grav- 
ity and coupled dark energy theories, dark energy evolves on 
cosmological time scales only when the scalar field leads to 
a long range interaction which could violate various gravita- 
tional bounds. To avoid this problem, screening mechanisms 
have been designed to dynamically screen the scalar-mediated 
fifth force in dense or high-curvature environments. 

Screening the effects of a scalar interaction in the presence 
of matter can be realised in the following ways. Let 0o be the 
environment dependent background configuration and let us 
expand the scalar Lagrangian to second order. 



SC = - 



{d5(t>f - 



where Z{(t)Q) is the normalisation of the scalar, to(0o) is 
the mass depending on the background value of the scalar 
field, and /?(0o) the coupling to the overdensity 6pm- The 
fifth force is screened if either /3((/)o) becomes small, or 
m(0o) or Z{(j>o) become large. Hence there are three known 
screening mechanisms to evade the local gravity constraints; 
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the Vainshtein mechanisrtiJQl in the Dvali-Gabadadze-Porrati 
(DO?) H and Gahleon models where Z{(j)o) is large 

enough to reduce the effective coupling /3(0Q)/Z^/^((/)n) be- 
low observational levels, the chameleon mechanism Q 01 
where the mass m(0o) is large enough to render the range 
of the scalar interaction smaller than distances probed exper- 
imentally and finally in dilaton [9] and symmetron ifiol ll ill 
theories the coupling ^{(pa) itself is smaller than observed. 
These mechanisms all utilise the nonlinearities of the effec- 
tive scalar Lagrangian to prevent the fifth-force from propa- 
gating freely. The nonlinearities for the Vainshtein case stem 
from the derivative self-couplings of the scalar degree of free- 
dom, while the chameleon, and the dilaton and symmetron use 
the non-linearities of the potential and the coupling to matter 
respectively. In a companion paper tl2il . we have presented 
a systematic study of generic dilaton and symmetron theo- 
ries in the nonlinear regime of structure formation, using A^- 
bodysimulations based on a unified parameterisation scheme 
iHlll.This paper will concentrate on chameleons. In partic- 
ular, we have generalised the original chameleon models and 
the parameter space of these new models is analysed making 
use of A^-body simulations. 

The structure formation in general chameleon models is dif- 
ferent from that in GR within the Compton wavelength of the 
scalar degree of freedom, which is the inverse of the effective 
mass of the scalar field m |[l4l - [l6ll . Indeed the density con- 
trast of matter perturbations increases anomalously there. This 
implies that the power spectrum differs from its GR counter- 
part. It turns out that the scale characterising modified grav- 
ity, i.e. the scalar mass now, is large enough to prevent sig- 
nificant effects on linear scales. The main consequences of 
modified gravity appear in the non-linear regime where nu- 
merical methods have to be used. The analysis of the struc- 
ture growth of the screened modified gravity models with no 
higher derivative terms in the Lagrangian, e.g. chameleons, is 
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rendered easier by the fact that these models can be fully pa- 
rameterised by two time-dependent functions, m(a) and /3(a), 
i.e. the mass of the scalar field and its coupling to matter as 
a function of the scale factor. This method works even on 
fully nonlinear scales |[l3i [Till where the screening effects 
on smaller scales, for instance for galaxy halos, can therefore 
be powerfully analysed. 

In this work, we study the nonlinear structure formation 
using A^-body simulations for the general chameleon theo- 
ries that we define here and which are parametrised by two 
m(a) and /3(a) functio ns 11311 . Technically, we use a variant 
of the ECOSMOG code 12311 . which is based on a public adap- 
tive mesh refinement (AMR) code RAMSES [|24l], to solve and 
evolve the A^-body system. 

Although chameleon theories, even in the nonlinear regime, 
have been studied extensively in the literature (e.g., ll25l - l40ll ). 
those studies are mostly for specific models in a very restricted 
parameter space. As an example, simulations for f{R) grav- 
ity have thus far only been done for the Hu-Sawicki model 
112611 which is equivalent to a chameleon theory with the cou- 
pling strength /3(a) fixed to l/\/6. Our study here, for the 
first time, allows /3(a) to have a time evolutior0.Our parame- 
terisation allows us to follow a more systematic approach to 
vary the different chameleon parameters and study the effects 
quantitatively. In particular, we find that the chameleon mech- 
anism is considerably more efficient than the dilaton and sym- 
metron mechanisms in restoring GR in high-density regions 
and at earlier times. Our results here show that linear perturba- 
tion theory fails almost whenever it predicts a deviation from 
ACDM, and point out the portion of the chameleon parame- 
ter space that is relevant to cosmology today and in the near 
future. 

The layout of this paper is as follows: in § HI] we review 
scalar-type theories and show how they can be parameterised 
simply; in §|III|we briefly describe the generalised chameleon 
model and the possible effects of varying each model parame- 
ter; the equations that will be used in the iV-body simulations 
are summarised in §|IVl various tests of our code are presented 
in §|V]and then the cosmological simulations used in this work 
are discussed in § IVII finally we summarise and conclude in 

§1211 

To make things clearer, throughout the paper we use the 
units h = c = 1 except where we use c explicitly. An overbar 
(subscript o) denotes the background (present-day) value of a 
quantity and subscript ^ means d/dip. k = SttGn = ^'^p^^ 
where Mpi is the reduced Planck mass and Gn is, Newton's 
constant, are used interchangeably for convenience. 



II. SCALAR-TENSOR THEORIES OF MODIFIED 
GRAVITY 

In this section we briefly describe the essential features 
of modified gravity theories with a scalar degree of freedom 
(dof) and how the effects of such a dof can be screened locally 
to restore general relativity (GR). More detailed descriptions 
can be found in our previous publications and here we keep 
the discussion short to make the paper self contained, and fa- 
miliar readers can skip this section. 



A. Scalar-tensor theories with screening 

The Einstein-Hilbert action for the scalar field (p in a 
generic scalar-tensor theory has the following form in the Ein- 
stein frame. 



S = 



Ml 




1. 



(1) 



in which g is the determinant of the metric tensor g^^ and R 

is the Ricci scalar. We label the ith matter field by -ipm . The 
quantities g^j, and g denote respectively the metric tensor in 
the Jordan frame and its determinant, and they are connected 
to g^i, and g via the following conformal transformation. 



5m'- = A'^{(p)g^^, g = A^{ip)g. 



(2) 



In the Einstein frame, the equation of motion (EOM) of the 
scalar field has an extra term because here tp explicitly couples 
to matter, and we get 



dip 



(3) 



in which T = —p + 3P is the trace of the energy momentum 
tensor T^", p, P are the energy density and pressure of matter, 
□ = V^V^j and the coupling strength between p and matter 
is given by fi{p) = M-pydlnA/dp. 

Eq. (O is equivalent to that of a normal quintessence field, 
with the bare scalar field potential replaced by a new effective 
potential 



Feff(^) ^ V{p) - {A{p) - 1)T. 



(4) 



In the simplest cases, V^s has a global minimum in the cosmo- 
logical background dominated by dust matter for which P,„ = 
and T = —pm- The value of the scalar field at the minimum 
depends on the actual value of p^, i.e., <^mm = Vm\n{Pm)- 
The mass of the scalar field at pn\\n, which is defined by 



' For other schemes to parameterise modified gravity see llT7H22r . Note how- 
ever that those schemes are mostly limited to the linear perturbation regime 
of modified gravity, while the parameterisation here is designed to account 
for nonlinearities. 

^ As we will see below, there are further subtle differences between the time 
evolutions of m(a) in our chameleon models and the model of t26il . 



m 



dp"^ 



(5) 



must be positive because an imaginary m can lead to violently 
unstable evolution of the perturbation of the scalar field. 
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When matter is described by dust fluid (with radiation neg- 
Ugible) so that there is no anisotropic stress, the Une element 
in the weak-field limit can be expressed as 



ds^ = -(1 + 2$)d<2 + (1 - 2$)dx2 



(6) 



where $ is the gravitational potential. This reduces to the 
modified geodesic equation for matter particles 



dV 



VM$ + ln^(^) 



(7) 



We can understand Eq. (|2| as the motion of a massive particle 
in an effective gravitational potential 



^off = $ + \iiA{ip), 



(8) 



and this is why the theory is considered as a modified gravity 
theory. 

As an example, let us consider a point mass M embedded 
in a homogeneous background density as the source of grav- 
ity. The effective gravitational potential could be obtained by 
solving the scalar EOM lfl2ll . as 



'I'off = - 



1 + 2/3((^)2e-™('^)'' 



GnM 



(9) 



The second term in the brackets represents a Yukawa-type de- 
viation from Newtonian gravity (the fifth force), and this devi- 
ation can be of order unity if mr < 1 and /3 ^ 0{1). However, 
because both /3 and m are functions of the field itself and thus 
depend on local matter density, it is possible that near massive 
bodies nonlinear effects make P{(p) ^ 1 or m^^ <C r. In 
such cases, the modification of gravity is strongly suppressed, 
which helps to evade local constraints on the fifth force. Be- 
cause the suppression of modified gravity here depends on the 
massive body itself, we call this ie//-screening. 

Self-screening is not the only mechanism to suppress the 
fifth force in modified gravity theories. Indeed, this suppres- 
sion often also very strongly depends on the environment of 
the body. In the case of chameleon theories, as shown in 13, 
the fifth force is effectively screened provided that the New- 
tonian potential <I>Ar at the edge of a massive body follows the 
relation 



Wo 



(10) 



where ipc is the minimum of Fcff inside the body and Lpoo , Poo 
are the minimum of Vcs and the coupling strength far away. In 
general, | (^c | ^ Woo], and <l>Ar in Eq. ( flOb determines the self- 
screening due to the massive body while ip^o (via also /3oo) 
characterises the environmental-scieemng. Note that although 
/3 is often chosen to be constant in chameleon theories, this 
does not necessarily have to be the case. 



B. Tomography 

As we shall see shortly, a rough estimate of the local con- 
straints on the fifth force indicates that rn^ ^ around the 



global minimum of Vcs- The scalar fiel d dy namically tracks 
</Jmin, around which it oscillates rapidly |[l4ll . and the time av- 
erage {Vcff (Vmin)) then acts as a very slowly-varying cosmo- 
logical constant. In this simplified case, we can determine the 
cosmic evolution of the scalar field in terms of to (a) and (3 {a) 
in background, as iflsilTill 



ip(a) 



il/pi am? {a) 



Pm{a)da + (pc 



(11) 



where we have assumed A{ip) = 1, as required by the strin- 
gent experimental constraints on the time variation of fermion 
masses, which is proportional to A. ipc is the scalar field value 
at the initial time Oini, when the average matter density in the 
Universe is of the same order as that in typical test bodies in 
laboratories today. Similarly, we have 



V{a) = Vo 



^. am? (a) 



Pm(a)da, 



(12) 



where Vq ~ V{a~ aini). 

Given V{a) and if{a), it is straightforward to derive V as 
a function of ip, V{ip). Similarly, l3{Lp) can be reconstructed 
easily from f3{a) and if{a). As a result, the full nonlinear dy- 
namics of the theory can be reconstructed elegantly using the 
background evolutions of to and /3. This 'tomography' ifisll 
has turned out to be very useful as a generic parameterisation 
of modified gravity theories and the systematic simulations to 
study their cosmological implications lfl2ll . 

We can then express the screening condition, Eq. (flOl i. as 



/9(a) 
am^ia) 



Prn{a)da < /3out A/pi's AT, 



(13) 



in which for simplicity we have considered constant matter 
densities Pout.in outside and inside the dense body, and out 
is defined by p™(ai„^out) = An^out and ^out = (3{a = flout)- 
One can use the fact that the Milky Way must be screeneqj 
to make a rough estimate about the screening condition. The 
averaged matter density inside the Milky Way is ^ 10^ times 
that of the cosmic mean, which implies that Oin ~ lO^'^; its 
Newtonian potential at its surface is ~ 10~^. On the other 
hand, approximately the environmental matter density for the 
Milky Way can be taken as close to the cosmic mearQ which 
gives us flout ~ 1- Using these numbers, Eq. (fTSl l implies that 
mo I Hq > 10'^. A similar bound can be deduced from the tim- 
ing of binary svste ms [|4lll and the distance indicators for stars 
in dwarf galaxies ll42ll . Hence, for a given modified gravity 
model to be screened locally, the fifth force can only act on 



^ It happens that the surface Newtonian potential is roughly the same for the 
Sun and the Galaxy, both ~ 0(10""). So if the Milky Way is not screened 
to provide environmental screening for the Sun, then the latter will not be 
self-screened either 

Clearly, this is only a simplified assumption, because the Milky Way lives 
in local high-density regions rather than the cosmological background. But 
here the purpose is only to roughly estimate the possible constraints coming 
from the Galaxy. 
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scales of and below a few megaparsecs in a cosmological set- 
ting. Because m itself is dimensional, in the rest of the paper 
we shall use the dimensionless quantity 



niQ 



(14) 



to parameterise modified gravity theories. ^ is proportional to 
the range of the fifth force, 



A = 2998^ h-^Mpc. 



(15) 



Even in GR, length scales of order megaparsec are already in 
the nonlinear regime and cannot be accurately described with 
linear perturbation theory. The nonlinearity in the equations 
for modified gravity only makes this situation even worse, and 
previous experiences ||40ll show that linear perturbation theory 
can be misleading whenever it predicts a deviation from GR. 
This has motivated us to analyse the large-scale structure for- 
mation in the chameleon theory more reliably, using A^-body 
simulations. 



III. GENERALISED CHAMELEON THEORIES 



A. Chameleon theory and its generalisation 



In the original chameleon theory proposed in |l7l|8|], the cou- 
pling function and the scalar field bare potential take the fol- 
lowing forms respectively: 



Vo 



A/pi 



(16) 
(17) 



Here /3o > is a dimensionless model parameter and Vq is a 
parameter with mass dimension four The chameleon screen- 
ing mechanism is graphically illustrated in Fig. [T] In high 
matter-density regions the contribution from the matter cou- 
pling to Vcff{(p) is large and the chameleon field (p is trapped 
in the small-field regime (i.e., — > 0) such that the fifth force, 
proportional to Vip, is very wealQ in low matter-density re- 
gions, in contrast, ip is big and so is Wip, resulting in a cos- 
mologically interesting fifth forc^ The essential features of 
the original chameleon theory include an exponential coupling 
function A{(p) and a runaway potential. 



^ As an example, in theories with a strong chameleon effect, the scalar field 
has very small value even in the background and under-dense regions, say 
i/p/Mpi £ [0, 10~*]. In this case, the variation of (p from the inside to the 
outside of a massive dark matter halo is at most ~ O(10~*), while the 
variation of the Newtonian potential is typically 0(10"^) or even larger, 
which means the fifth force is much weaker than standard gravity. Indeed, 
the smallness of is a generic consequence of the chameleon effect. 

* One can also understand the suppression of the fifth force in high matter- 
density regions as a result of the locally very heavy scalar field mass, which 
characterises the length scale the scalar degree of freedom could propagate 
without being severely suppressed. 



As discussed above, a coupled scalar field, if heavy enough 
(namely m(a) ^ H), can have its time evolution fully speci- 
fied by m{a) and /3(a), both of which are determined as func- 
tions of a by the background cosmology. For the chameleon 
theory listed in ifisll , it has been shown that 



m{a) 
/3(a) 



moa 
/3o. 



(18) 
(19) 



where r > 0. 

As a straightforward generalisation of the chameleon idea, 
in this paper we shall consider a power-law form of both m(a) 
(as in Eq. ( fTST i) and /3(a): 



/3{a) = /3oa 



(20) 



where s is a new model parameter to describe the generalised 
chameleon theory. Using the tomographic mapping discussed 
above, we find that 



il/pi 
A/pi 



1 



2r-s-3 



,2r-s-31 



2r- s-3 



where we have used a subscript i to denote the value of a quan- 
tity at the initial time a^, and ^ — Ho/mo as defined above. 
As we take the limit ai — > 0, the chameleon field is driven to 
(f ^ Q and the above equation reduces to 



Mpi 2r 



2 2r- 



(21) 



In order to study the nonlinear evolution of (p, we have to 
know Vip{p), where a subscript ^ denotes derivative with re- 
spect to p, which governs the dynamics of the scalar field (see 
the A^-body equations below). For this we find 



d{KV{a)) da 
da dp 
—3i}m/3oH^a 



s-3 



(22) 



— — 30,„/3oi/o 



9i^m/3o 


3+s 
2r-s-3 




2r - s - 3 







2r-s-3 

(23) 



in which Eq. ( l22b can be used in background cosmology and 
linear perturbation analysis to replace V^p. As /3o, ilm and 
are all positive, to make sure that the quantities in Eq. (1231 ) are 
well defined we will require 2r — s — 3 > and > in what 
follow^]. We also require that r > 2, since if cx a during 
the matter-dominated era and oc a~'^ in the radiation era, 
so that one may have > at early times if r < 2. 

With V^{p), one could easily integrate to obtain V{p) an- 
alytically and we have 



K,V{p) = kVq 



2r - 2s - 6 



2r 



3 p 



9nM^ Afpi 



(24) 



^ Otherwise the terms in the brackets could be negative, making the power- 
law function ill defined. 
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FIG. 1. Illustration of how the chameleon mechanism works. The dashed, dotted and solid curves are respectively the bare potential V{ip) of 
the chameleon field, the coupling function and the total effective potential Veff (v)- Left Panel: in high matter-density regions the minimum of 
Vc{i{(p) is very close to = and V</3 is small so that the fifth force is strongly suppressed. Right Panel: in low matter-density regions ip and 
therefore V</5 can be big, and so a nonzero fifth force takes effect in structure formation. 



for r — s 7^ 3 and 



27, 



KViif ) = kVo - -nt^l3^eH^ log 



PI 



(25) 



for r — s = 3. The perturbation of the dark energy density, 
V{ifi) — V{(p), appears in the source for the Poisson equation 
(see below), but it is generally very small and can be safely 
neglected. 

Given ip{a) and /3(a), it is straightforward to find (3{ip). For 
our parameterisation using r, s, we find that 



A{^) = 1 



/3o 



2r-s-: 



2r- 



2r - 2s - 3 



and 



Mpi 



(26) 



2r- 



3 (p 



(27) 



As a result, both and (3 are power-law functions of (f. 



B. Effects of varying chameleon parameters 

As shown above, our generalised chameleon theory is spec- 
ified by four model parameters, namely, f3o, r, s and ^. The ef- 
fect of varying these four parameters on the structure forma- 
tion can be understood without solving the system explicitly. 

The parameter /3o, which is the coupling strength at z = 0, 
controls the overall amplitude of the coupling throughout the 
entire evolution history. The larger /3o is, the stronger the fifth 



force is, thus the strong clustering of matter relative to that in 
ACDM (in which = 0). 

The parameter r, which is the power index of m{a), deter- 
mines the time evolution of the effective mass of the scalar 
field without changing niQ, which is the mass at z = 0. The 
smaller r is (?- > 0), the lighter the scalar field is at 2 > 0, and 
so the longer the range of the fifth force is. Due to the tomog- 
raphy mapping, this also implies that the scalar field is less 
heavy in high-density regions, leading to a weaker chameleon 
screening and a stronger clustering of matter Recall from the 
above that we restrict ourselves to r > 2. 

The parameter s, which is the power index of /3(a), deter- 
mines the time evolution of the coupling function. The more 
negative s is, the weaker the coupling between matter and the 
scalar field at z > becomes; because of the tomography re- 
lation, this also means a stronger suppression of the fifth force 
in high-density regions, and therefore weaker matter cluster- 
ing. Note that here we restrict ourselves to s < to avoid the 
anfi-chameleon effect: this can be seen by looking at Eq. dZTl l. 
which shows that the coupling is stronger in high density re- 
gions, or Eq. (|20] i. which shows that the coupling is stronger at 
earlier times. The situation is worse if ?■ — 3/2 < s < 2r — 3, 
which implies that A{ip) decreases with ip by Eq. (l26b , and 
there is no longer any minimum for VcB- 

The parameter ^, which is simply Ho/mo, essentially sets 
the effective mass m of the scalar field (and thus the effective 
range of the fifth force) at z = 0. In all the chameleon simula- 
tions we study in this work, ^ ^ 1. The larger ^ is, the lighter 
the scalar field is and the stronger the fifth force becomes. 

In what follows, we will find that the A^-body simulations 
confirm this analysis and also quantify these effects. 
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IV. THE iV-BODY EQUATIONS 

This section serves to introduce the A^-body Poisson and 
chameleon equations for the sake of completeness. For this we 
list the equations to be solved and describe the code units used 
in our simulations, both of which can be found in 



A. Simplified field equations 

The relevant equations which determine the dynamics of 
the chameleon and gravity fields are 



dV 

dt2 



AttG {pm - Pni) , 



(28) 

A^i^)p,n,(29) 



dr 



(30) 



where we work in the quasi-static limit by dropping all terms 
involving time derivatives. 

The validity of the quasi-static approximation was tested 
explicitly in ll28ll . which compared the time and spatial deriva- 
tives and found that the former is indeed negligible. Note that, 
rigorously speaking, |i28|] only tested that ^ (cp) is negligible, 
where (ip) is the scalar field value averaged over the quick os- 
cillations, rather than ip itself, which can be as large as | V(^| 
due to the oscillations. The oscillations themselves, however, 
largely cancel out and it is the averaged effect that we observe 
- in this sense we believe that the test of ll28ll is accurate. We 
have checked, using our linear perturbation code, that the ef- 
fects on cosmological observables (such as as) differ by less 
than ^ 0.1% in the two cases where we respectively follow 
the oscillations accurately and average over them l!l4ll . 

It is tem ptin g to try to solve the full time-dependent scalar 
field EOM [43i] in modified gravity simulations, but notice that 
to follow the time evolution one has to resolve the oscilla- 
tions very well. It does not seem so difficult at the background 
level, where mo/Ho ^ lO'^, meaning that to accurately re- 
solve the oscillations one needs a factor of 0(10) xe'(103) - 
0(10"*) — 0(10^) coarse time steps. However, evenin a mildly 
high-density region one could have miocai/^^o ~ 10®, re- 
quiring 0(10^) — 0(10^) time steps to accurately follow the 
time evolution during the course of an iV-body simulation. 
For comparison, the simulations in this paper uses a few hun- 
dred coarse time-steps so fully solving the EOM represents a 
huge increase in the computational cost of a simulation. Us- 
ing fewer time steps would mean that some sort of average 
has been done implicitly, in the same sense as it is done in the 
quasi-static approximation. 

A full treatment of this issue is beyond the scope of this 
work, but rigorous test of the quasi-static approximation for 
modified gravity theories is something we plan to pursue in 
the future. 

B. Code units 

The code units are based on (but not exactly) the superco- 
moving coordinates of lliill . They can be summarised as fol- 



lows (tilded quantities are expressed in the code unit): 



X 



Ha^, c 



pa 



Pcflr. 



av 



where x is the physical coordinate, t is the physical time, c is 
the speed of light, pc the critical density at present, ftm today's 
fractional energy density for matter, v the particle velocity and 
<i> the Newtonian potential. Besides, B is the simulation box 
size in unit of h^^Mpc. The average matter density is p = 1 
in the code unit. Note that all these code quantities are dimen- 
sionless. 



C. The discrete equations 

In cosmological simulations, the chameleon field p is gen- 
erally extremely small {i.e., p/Mp\ <C 1) and must be positive 
to make the logarithmic in Eq. ( |23] | well defined. To prevent 
the numerical value of p from becoming negative during the 
computation and therefore causing divergence problems, we 
follow ll28l[3ll [33I1 to define a new variable u = \og{p / Mpy) 
instead of using p itself. Throughout the cosmic evolution and 
from one spatial position to another, p can change by several 
orders of magnitude, but |u| remains 0(1 ~' 10), making the 
numerical problems easier to avoid when using u. 

Using the quantities defined above, the Poisson equation, 
Eq. (|28l) . can be written as 



(31) 



and, after some manipulation, the chameleon equation of mo- 
tion Eq. ( |29l ) reads 



V"e 



2u 



^VtmPap 



e 

^ . 

L <^ 



'1~S 



(32) 



Before being implemented into the iV-body code, the above 
equations must be discretised. For the Poisson equation this is 
straightforward and we have 



1 

7? 



+ ^ij^k-l - 6$ij,fc] = -ilrnCL {Pi,j,k - 1) (33) 



where 'tij^k denotes the value of $ in the {i, j, k)-th cell of the 
simulation grid. The discretised nonlinear chameleon equa- 
tion can be obtained in a similar way though it involves longer 
derivation. 



L''(w,,,-fc) =0, 
with the operator L'^{ui j k) defined as 



(34) 
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(35) 



Here h = de^/du = e", 

and /i is the cell size in the simulation grid0. 



In our simulation, Eq. (|34| ) will be solved using the Newton 
Gauss-Seidel relaxation method, described as 



h,ncw h,o\d V ij^k J 



where 



(36) 



dL'' {Ui fe) £\ r . 1 



+ bi-i,j.k + bij+i^k + bij-i.k + bij,k+i + bij^k-i + 66ij-,fc] 



2r~Su 



,2r-3 



2r - 3 u. 



(37) 



TABLE I. The parameter values for the seven models used in the 
code test. 



parameter 


/3o 


r 


s 




model a 


0.5 


3.0 


0.0 


0.001 


model b 


1.0 


3.0 


0.0 


0.001 


model c 


0.5 


4.0 


0.0 


0.001 


model d 


0.5 


3.0 


-1.0 


0.001 


model el 


0.5 


3.0 


0.0 


0.0005 


model e2 


0.5 


3.0 


0.0 


0.002 


model e3 


0.5 


3.0 


0.0 


0.005 



A. Homogeneous matter density field 

In a homogeneous matter density field, the chameleon field 
ip must take a constant value given by 



^ = ^ ^ S lM^c?''-^-^ (38) 
2r — s — 6 



V. CODE TESTS 

To make sure that our code works properly, we performed a 
number of code tests, which are described in this section. We 
tested the code for 7 models by varying the 4 parameters for 
the generalised chameleon model, namely (3o,r,s and ^, and 
these are summarised in table |V] Throughout this section we 
adopt the unit Mpi = 1. 



Note that h is also used in this paper to denote _ffo/(100 km/s/Mpc), 
but there should be no confusion since it is easy to understand its actual 
meaning based on the context. 



Our first test, therefore, is to fix the matter density field on the 
simulation grid, make a random initial guess about the values 
of u, let the Newton Gauss-Seidel relaxation iterate for a few 
steps, and see if u approaches \ogLp [ip given in the above 
equation) in all grid cells. 

We did this test for 3 out of the 5 models described in Ta- 
ble |V] at a = 1.0, as shown in Fig.|2l in which we plotted the 
values of u in all cells in the x-direction, both before (open 
symbols) and after (filled symbols) the relaxation. We could 
see clearly a good agreement between the numerical solutions 
(filled symbols) and analytic results (the horizontal lines). We 
also did the test at several values of a < 1.0 and found similar 
agreements, but these are now shown here for clarity. 
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FIG. 2. (Colour online) Homogeneous matter field test. Shown are 
the values of the scalar field (p in all cells along the a::-direction, with 
y = 2 = 0. To make the plot clearer, only results for models a and b 
are shown. Open symbols represent the initial guess and filled sym- 
bols of the same shape and colour represent the numerical solutions. 
Please note that, instead of itself, shown here is it = log(<^), which 
is what the code outputs directly. The thick solid lines are the analyt- 
ical solutions specified in Eq. ([38} for the two models, and they have 
the same colour as the corresponding numerical solutions. 



B. Point mass 

Our second test makes use of the simplest spherical ly sy m- 
metric density field, a point mass at the origin ll23ll28[l5Tll . In 
which case there is an exact analytic solution to cp (or equiva- 
lently u) under certain simplifications. 

The said density configuration can be constructed ll28ll as 



10-* 
-10- 



{N^ - 1) 



i ^ j = k = 0; 
otherwise. 



(39) 



in which we have defined Sij,k = Pij.k — 1- The analytical 
solution can be obtained by solving the equation 



(40) 



away from the origin (where the point mass is), in which the 
mass of the chameleon dof, S^p = — (p, is given by = 
£,^Hq, and by doing that we found 



Sip (X — exp(— mr), 
r 



(41) 



in which r is the distance to the origin. 

In this test, the simulation box we chose has a length of 
250/i^^Mpc and 256 cells on each side, and we tested all the 
5 models at a = 1.0. Fig.[3]compares the numerical solutions 
(symbols) of Sip in the x-direction to the analytical predic- 
tions (solid curves) given in Eq. ( |4TI ). and it is clear that they 
agree very well for all tested models. We stress that the dis- 
crepancies on large and small values of r are not indications 



FIG. 3. (Colour online) Point mass tests. Shown here as filled sym- 
bols are the values of Sip = ip — ip away from a point particle con- 
structed as described in Eq. ( 139) , along the x-direction. To make the 
plot less busy, we only show the results for 4 out of the 7 test models 
listed in Table IVl (more details in the legends). We use solid curves 
for the analytical solutions of Eq. ( I41t . which are good approxima- 
tions if the distance to the point mass is not too small. 



of the code's failure - the former is because the magnitude of 
Sip is already at the level of discretisation error (which limits 
the code's ability to make the solution more accurate by fur- 
ther relaxation iterations, and this can be seen from the fact 
that the discrepancy happens at the same value of Sip for all 
models), and the latter is because of the fact that in deriving 
Eq. ( |40b one artificially Hnearises a nonlinear equation Ii28il . 



C. Sine density field 

The next two tests make use of one-dimensional density 
configurations, which are obtained by eliminating the y- and 
z-dependences of the density field. Starting from a given ID 
solution to ip, we substituted it into the chameleon EOM to 
find the desired density field, and then used this density field 
in the numerical code to solve for ip and compare with that 
original input. 

The first such test uses a sine density field as first introduced 
in 12811. which in our code units can be written as 



f [2 - sin(27ra;)] 



siii(27ra;) [2 — sin(27rx)]- 



(42) 



in which x is rescaled by B and so a; € [0, 1]. The correspond- 
ing ip field which is associated wth this density field is 



^p{x) = (ys [2 — sin(27r2;)] . 



(43) 



We did this test for 6 models listed in Table IVl at a = 1.0 
and the results are shown in Fig. |4l As expected, there is a 
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FIG. 4. (Colour online) Tests with sine density fields, for 6 out of the 
7 test models (a, b, c, d, el, e2) at a — 1.0 (see the legends). Shown 
here are the numerical solutions (filled symbols) of ip along the x- 
direction in the sine-type density field described in Eq. (142b . and the 
analytical results of Eq. ( I43t (solid curves of the same colour as the 
symbols). For this test the simulation box size is 250/i"^Mpc and 
X is rescaled to make x/B £ [0, 1]. The simulation mesh has 256"^ 
cells. 
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FIG. 5. (Colour online) Tests with Gaussian density fields, for 6 out 
of the 7 test models (a, b, c, d, el, e2) at a = 1.0 (see the legends for 
details). Shown here are the numerical solutions (filled symbols) of 
ip along the a;-direction in the Gaussian-type density field described 
in Eq. i45\ . and the analytical results of Eq. ( I44t (solid curves of the 
same colour as the symbols). For this test the simulation box size is 
250/1" ^Mpc and X is rescaled to make x/B £ [0, 1]. The simulation 
mesh has 256^ cells. 



good agreement between the numerical solutions (open sym- 
bols) and the analytical results (filled symbols). 



D. Gaussian density field 



Fig. |5] shows the test results for 6 out of the 7 models sum- 
marised in Table |V] at a = 1.0, and again the numerical solu- 
tions (symbols) match the analytical solutions (soUd curves) 
of Eq. ( l44b ver accurately. 



The second test that makes us of a ID density field assumes 
a Gaussian-type solution to (p, given by 



1 — a cxp 



(a; -0.5)^ 



(44) 



where W, a are numerical constants which are used to specify 
the width and height of the Gaussian function. As before, x is 
scaled by the boxsize so that x G [0,1]. Note that the Gaussian 
function in ip{x)peaks at x = 0.5 while at a; or .t 1 we 
have If ^ (p. Also, a 1 makes \ip\ very small at x = 0.5. 

The density field which is associated to the above solution 
to (p{x) is 



p{x) 



{x~o.5y 



-L — ^ 



1 — a exp 



(x-O.S)^ 



1 



ae 



(45) 



Notice that such a density field is not exactly periodic at the 
edges of the simulation box, but given that W is small enough, 
p — > at the box edges and periodic boundary conditions are 
approximately satisfied. 



E. Multilevels 

One of the most important features of the ECOSMOG code 
is that it enables adaptive mesh refinements, and to make sure 
that this part of the code also works correctly we need to test 
it on the refinements. This is the task of this subsection. 

The Gaussian test described above provides a good starting 
point for the multilevel test here, because the density contrast 
at a; = 0.5 could be made very large by choosing appropri- 
ate values for a, which triggers refinements of the simulation 
meshes. Indeed, when a 1, the fast change of density field 
close to a; = 0.5 makes refinements essential to guarantee the 
high precision. For simplicity, in the test here we only refined 
the grid once, making this a 'two-level problem', with level 
8 (9) representing the coarse mesh (refinement), where 'level 
8' means the mesh has 2* = 256 cells in each dimension. On 
both levels we used Eq. ( |45] ) to set the density values in cells, 
and for level 9 we set the boundary conditions for ip by in- 
terpolating the corresponding values in the coarse cells which 
cover the refinement boundary (more details can be found in 
Q). 

Fig.|6]shows the test results for model a only and for 3 dif- 
ferent values of a (0.999, 0.9999 and 0.99999 from top to 



10 



1E-10 



1E-11 



_ 

: 

- 


1 1 1 1 1 1 1 1 1 1 


1 1 1 1 : 

. 

- 


r 


" a=0. 999, level 3 \ t 






■ 0=0.999. level 9 \i 






a=0.9999, level 8 \ i 






0=0.9999, level 9 






0=0.99999, level 8 






• 0=0.99999, level 9 " 




1,1, 



0.42 0.44 0.46 0.48 0.50 0.52 0.54 0.56 0.58 



x/e 



model name /3o r s ^ realisations 



ACDM _ _ _ _ 5 

baseline 0.50 3.0 0.001 5 

Al 0.25 3.0 0.001 5 

A2 0.75 3.0 0.001 5 

Bl 0.50 2.0 0.001 5 

B2 0.50 2.5 0.001 5 

B3 0.50 3.5 0.001 5 

CI 0.50 3.0 -0.25 0.001 5 

C2 0.50 3.0 -0.5 0.001 5 

C3 0.50 3.0 -1 0.001 5 

Dl 0.50 3.0 0.0005 5 

D2 0.50 3.0 0.0015 5 

D3 0.50 3.0 0.002 5 



TABLE II. The parameter values for the 65 cosmological simula- 
tions we have performed for this study. Note that '-' means that the 
parameters are unused for the ACDM case. 



FIG. 6. (Colour online) Multilevel tests for model a at a = 1.0 and 
three different values of a: 0.999 (red), 0.9999 (green) and 0.99999 
(blue) from top to bottom. The open and filled symbols, with same 
shapes and colours, represent respectively the solutions to Lp along 
the a;-direction on level 8 (the domain grid) and level 9 (refinement), 
and the corresponding analytical solutions of Eq. (I44t are shown us- 
ing solid curves of the same colours. For this test the simulation box 
size is 250/i^^Mpc and x is rescaled to make xjB £ [0, 1] (code 
unit). The level-8 simulation mesh has 256"^ cells. 



bottom). In each case, we represent the numerical solutions on 
levels 8 and 9 by open and filled symbols of the same shape 
and colour, and the analytical results Eq. (|44] | by solid curves 
of the same colour Not surprisingly, the numerical and ana- 
lytical solutions agree very well; so do the numerical solutions 
on the two different levels. 

These tests make us confident about the reliability of our 
code, and about the simulations we describe below. 



VI. COSMOLOGICAL SIMULATIONS 
A. Simulation details 

This section is the core of this paper, and it shows the results 
of the cosmological simulations of our generalised chameleon 
models. We simulated a total 13 models with different values 
of /3o, r, s and ^, including the special case of ACDM which 
corresponds to /3o = 0.0, as summarised in Table IvTaI For 
each of the models we have 5 realisations to be averaged over 
to make the physical predictions more statistically meaning- 
ful, and all these 5 realisations have the same physical and 
simulation parameters, with the only difference being in their 
initial conditions, which are generated by MPGRAFIC ||45l] at 
an initial redshift Zi = 49.0 using different seeds of random 
numbers. 

The cosmic expansion rate in all our simulated chameleon 
models is very close to that of the standard ACDM paradigm 



lfl4ll . and this is determined by the WMAP7 ||46t] cosmological 
parameters. In particular, 

{h, n,n, f^A, Us, as} = {0.71, 0.267, 0.733, 0.963, 0.801}. 

Our simulation box is 128/i^^Mpc in each dimension, and the 
domain gric0has 256'^ cubic cells. Any cell is refined and split 
into 8 son cells when the number of particles inside it exceeds 
9.0, and in our simulations the finest refinement level has 2^^ 
cells on each side assuming that it covers the whole box. We 
use Np — 256'^ dark matter particles in the simulations. 

In the rest of this subsection, we will focus on the effects 
of changing each model parameter on the major cosmological 
observables such as the matter power spectrum and halo mass 
function. More specifically, we will analyse the results of our 
numerical simulations according to the following: 

1. How the amplitude of the coupling strength today, /3o, 
affects the results: models Al and A2; 

2. How the power of the scalar field mass, r, affects the 
results: models BI, B2 and B3; 

3. How the power of the coupling strength, s, affects the 
results: models CI, C2 and C3; 

4. How the range of the fifth force, ^, influences the re- 
sults: model Dl, D2 and D3. 

To see more clearly the effect of varying these four parame- 
ters, we have also simulated a baseline model {/3o, s, = 
{0.50, 3.0, 0, 0.001}, to which all other models are compared. 

B. Nonlinear matter power spectra 

The most direct way to see the effect of modified gravity on 
the clustering of matter is to look at the matter power spectrum 



' In AMR codes such as RAMSES and ECOSMOG, the domain grid is the 
unifomi (regular) grid which covers the whole simulation domain. 
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FIG. 7. The fractional difference in matter power spectra of various chameleon models (different models are illustrated in the legend) with 
respective to that of the ACDM model at z = 0. The curves with error bars show the simulation result, while the curves without error bars 
stand for the linear theory prediction. In each panel, the curves with the same color and line style represent the same chameleon model. 



P{k). For this we have measured the P{k) for our generalised 
chameleon theories and the ACDM paradigm using POWMES 
ll47ll . and calculated the relative difference AP/Pqr - The re- 
sults are shown in Figs.|2]and[8] 

In Figs. |2] and [8] we can see that both the linear perturba- 
tion results (the smooth curves) and the simulation predictions 
(symbols) follow the trend as we have expected (see § 1111 BI) . 
The linear perturbation prediction significantly overestimates 
the relative growth with respect to that in ACDM model in 
all cases, similar to what we found in the dilaton, symmetron 
and f{R) gravity simulations lfT2[|40ll . In particular, we notice 
from these figures that, linear perturbation theory fails when- 
ever it predicts a deviation from ACDM, and this can happen 
on scales as large as /c ~ O.OSMpc^^. This result casts strong 
doubts on all the efforts which have been made to constrain 
chameleon-type theories using linear theory predictions, and 
shows once again the crucial role played by nonlinear simula- 
tions. 

It may seem to be surprising that linear theory breaks down 
on large scales which can be well described by it in standard 
cosmology, and the reason is that the chameleon theory itself 
is nonlinear, and this nonlinearity is in addition to the usual 
nonlinearity in real matter distributions. Consequently, in the 
fifth-force calculation different Fourier modes of the density 



field strongly couple and the fifth force for large-scale modes 
depends on the matter perturbation on smaller scales. In linear 
perturbation theory, such mode coupling has been suppressed. 
Another way to understand the point is the following: in linear 
theory, the fifth force (is assumed to) depend only on the back- 
ground matter density, while in nonlinear simulations it actu- 
ally depends on the matter density inside overdensities which 
is generally higher than the background density and therefore 
makes it more suppressed due to the chameleon mechanism. 

Note that the agreement between linear perturbation theory 
and A^-body simulation results is up to smaller scales at 2: = 1 
than at z = 0, but this is most likely because both approaches 
predict smaller deviations from ACDM at earlier times when 
the matter density is higher overall, rather than because linear 
theory works better at higher redshifts when density perturba- 
tions are small. Indeed, a direct comparison between Figs. [T] 
and [8] confirms that the (nonlinear) chameleon effect is much 
stronger at early times. 

The upper left panel of Fig. |7] shows the effect of varying 
f3o while all the other parameters are fixed to their baseline 
values (c.f. Table I VI Ab . As shown, AP/Pqr increases when 
f3o rises. Specifically, we increase and decrease /3o around 
0.5 (which is the value of the baseline model) by 50% in 
models Al and A2 respectively, and find strong variations in 
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FIG. 8. The same as Fig|7] but al z — 1. 



the linear theory predictions of AP/Pqr- The simulation re- 
sult of AP/Pqr, however, is smaller than ^ 2% down to 
k = even at z = 0. This small deviation is be- 

yond the precision of all current cosmological probes. Recall 
that /3o here is chosen to have the same value as that in the 
dilaton simulations of lfl2ll . where AP/ Pqr can be more than 
~ 30 — 40% - this shows clearly that the chameleon screening 
is much more efficient in restoring GR in dense regions. 

The upper right panel of Fig. [T] shows the effects of vary- 
ing r while other parameters are all fixed to the baseline val- 
ues. The result is again consistent with the analysis in § 1111 B I 
namely, AP/ Pgr grows as r drops because a smaller r means 
a less massive scalar field in the past or, thanks to the tomog- 
raphy mapping, in dense regions. For example, the chameleon 
screening in model Bl is less efficient than that in B3 at z > 0, 
making gravity relatively stronger in the former during most 
of the the evolution history, which is why the accumulated ef- 
fect on matter clustering is much more significant in B 1 . 

The lower left panel of Fig. [T] illustrates the effect of vary- 
ing s while other parameters are fixed to the baseline values. 
As expected, AP/ Pqr drops as s decreases, which is because 
a smaller coupling in the past or in dense regions necessarily 
means a weaker fifth force and therefore a decrease in the mat- 
ter clustering. As we mentioned above, to avoid the unwanted 
anti-chameleon effect we have to choose s < 0, which means 
that the baseline model, with s ~ 0, gives the largest possible 



deviation from ACDM, which is < 1% at fc = 0.1/iMpc^^ - 
this clearly implies that s is practically unconstrained except 
that s < 0. 

Finally, in the lower right panel of Fig. |2]we have shown 
the effect of varying ^ with all other parameters fixed. Since 
^ is inversely proportional to niQ, an increase in ^ results in a 
smaller scalar field mass throughout the evolution history and 
therefore more structures form due to the weaker suppression 
of the fifth force. This is exactly what we see in this panel. 

Overall, Figs.|2]and[8]indicate that observational data on the 
matter clustering at present and in the near future will hardly 
place any strong constraints on the chameleon-type modified 
gravity theories. One therefore has to look at other cosmolog- 
ical probes, such as the halo mass functions and void proper- 
ties, to detect any observable signatures of these theories. We 
will study the former in the next subsection and leave the latter 
to future work. 



1. Comparison witii f{R) gravity model 

Note that the models we study in this work generally have a 
much stronger chameleon effect compared to the f{R) models 
simulated in Isil l40ll . which are the Hu-Sawicki model ll26ll 
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with n = 1 and \ fRo\ = 10-^ 10-^ 10"'' respectivel>B 



From Eqs. (12, 13, 18) of 134||, it is straightforward to find 



m 





(<■■ 




Of ^ 


2\fm\ 1 




1 



(46) 



From this expression we can immediately learn two things. 
First, the ^ parameter in the Hu-Sawicki f{R) model is given 
by 



mo 



2|/. 



R0\ 



(47) 



Taking f2,„ = 0.25, f^A = 1 - r2,„ = 0.75 and Jro = -10"^ 
we have ^ w 0.78 x 10^"^. Second, m{a) is a power-law func- 
tion 



m(a) oc a 



-4.5 



(48) 



withr = — 4.5fora~^ > 3, while for 0(1) then m(a) 

stays almost a constant. In addition to these, it is well known 
that f{R) gravity is a special case of chameleon theories with 
/3o = 1/v^ands = 0. 

Judging form the values of ^, s and /3o, it may seem that the 
Hu-Sawicki model with fm = —10^^ should lead to smaller 
deviation from ACDM than the baseline model. It looks even 
more so if one considers that r — —4.5 < —3 for small a, and 
this seems to be inconsistent with the simulations. Note here, 
however, that r = —4.5 only happens for z ^ 1 when the 
fifth force is negligible anyway, and at z < 1 m stays around 
7710 so that the fifth force is indeed less suppressed than in the 
baseline model here. 



C. Dark matter halo mass functions 



while at z = 1, smaller halos with mass M > 10^^h~^ Mq 
can also be well screened in some, if not all, cases. 

The effects of varying the different chameleon parameters 
are generally the same as what we have expected or have seen 
in the plots of AP/Pgr- namely, Au/hq^ increases as f3o,s 
and ^ increases or r decreases. Different from the case of mat- 
ter power spectra, however, the mass functions in chameleon 
theories show larger deviations from that of ACDM, particu- 
larly in the low-mass end. 

A nontrivial feature in Fig. |9]is the turnover on An / ncR for 
models Bl, B2, Dl and D2. Without loss of generality, let us 
take model D2 as an example and com pare to the f{R) model 
with |/_Ro| = 10~^ (F6) simulated in Il34ll . In both cases, the 
largest halos in the simulation box are well-screened, both by 
themselves and by their environment (because large halos tend 
to be produced out of very dense regions). When the halo mass 
decreases, the self-screening becomes weaker and the halo has 
a higher probability of living in average, or even underdense, 
regions - the weakened screening means more matter cluster- 
ing and production of more halos. Of course, there is a limited 
supply of matter to be incorporated into halos, and when more 
large halos are formed there will be fewer small halos surviv- 
ing the mergers and accretions, that has caused the turn-over 
This is the same as what is found for the F6 model in ll34ll (see 
Fig. 1 1 therein) and also complies with the analytical results 

of iilll. 

The chameleon effect in the rest of our simulated models is 
too strong so that even low mass halos get screened to a certain 
extent, making the growing trend with mass at the low mass 
end disappear. This can be seen by looking at the D2 model in 
FigHni the turnover disappears simply because the chameleon 
is more efficient at higher redshifts. Also note that at z = 1 
the suppression of the fifth force is so strong that the deviation 
from ACDM almost vanishes for most models, which is the 
same as we have seen in the AP/Pqr plots above. 



We measured the dark matter halo mass functions from our 
simulations using the publicly available code AHF [48i], which 
is efficiently parallelised using MPI and OPENMP. We define 
the halo mass as the total mass contained in i?200- the radius 
at which the average matter density inside drops below 200 
times the critical density. For each model, we have calculated 
the binned relative difference in mass function with respect to 
that of the ACDM model (see IH] for details). 

In Figs.|9]and[T0]we show the ratios between the chameleon 
and ACDM mass functions from our simulations at z = and 
z = 1 respectively. From these figures it can be seen clearly 
that the fifth force leads to an overall enhancement of the for- 
mation of dark matter halos, and the effect is stronger on the 
low-mass end of the mass function. The maximum An/nQR 
is around 50% for the models we simulated. At z = 0, halos 
with mass M > 5 x 10^^h~^MQ are generally well screened. 



For more details of the models and the definitions of fno and n, see flS] or 
I34II4OII . Here we will quote the results rather than give a thorough review. 



VII. SUMMARY AND CONCLUSIONS 

To summarise, in this paper we have brought together two 
essential techniques for the systematic studies of the nonlin- 
ear structure formation in generic modified gravity theorie s of 
the chameleon type: a simple parameterisation scheme which 
covers all known chameleon theories using only four parame- 
ters and a modified version of the ECOSMOG code to run high- 
resolution simulations efficiently. This allows us, for the first 
time, to get an overall picture about the behaviour of general 
chameleon-type theories and the part of its parameter space 
which is relevant for cosmology. 

The powerful tomography mapping lfl3l [Till enables us to 
characterise the chameleon theory and its generalisations us- 
ing only a few parameters. In our case, there are two parame- 
ters describing the present value of the scalar field mass 
and its time evolution (r), and another two parameters de- 
scribing the current value of the coupling strength (/3o) and its 
time evolution (s). These 4 parameters cover most chameleon 
theories studied in the literature jljl . and also the cases with 
varying (field-dependent) coupling to nonlinear structure for- 
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FIG. 9. The fractional difference in halo mass function of various chameleon models (different models are illustrated in the legend) with 
respective to that of the ACDM model at 2; = 0. 



mation which have not been thoroughly investigated so far. 

Following the logic of |12], here we focus on the qualitative 
and quantitative behaviour of the generalised chameleon the- 
ory. We are interested not only in how varying the parameters 
changes the predictions of cosmological observables, but also 
in how large the changes could be such that we can decide 
which portion of the 4D parameter space would be of inter- 
est to cosmologists and therefore merits further (and more de- 
tailed) investigations in the future. As a by product, we want 
to compare the efficiencies of the different screening mecha- 
nisms that have been explored by theorists - the chameleon, 
dilaton and symmetron mechanisms. 

To this end, we have simulated a total of 12 models which 
form an extensive span in the parameter space. Starting from 
adefaultmodel with {/3o, r, s, ^} = {0.5, 3.0, 0.0, 0.001}, we 
let each of the 4 parameters vary and take a few different val- 
ues as summarised in Table IVI Al In this way, we can see 
clearly the effect of changing every parameter. 

The simulation results confirm our qualitative predictions 
based on simple physical arguments, namely the the fifth force 
(and therefore the clustering of matter) is stronger if one; 

1. increases /3o, which results in an overall increase in the 
coupling strength between matter and the scalar field; 

2. increases s, which makes the coupling strength reduce 



more slowly as the matter density increases; 

3. increases ^, which increases the range of the fifth force 
overall, or 

4. decreases r, which makes the fifth force less exponen- 
tially suppressed in high-density regions. 

There are a few noticeable features which can be seen from 
the nonlinear matter power spectrum predicted by our simula- 
tions. The first is that, as in the cases of dilaton 11211 and f{R) 
gravity | |40l ] models where the screening is strong, linear per- 
turbation theory fails for general chameleon theories wherever 
it predicts a deviation from ACDM. The scale at which linear 
theory breaks down can be as large as fc ~ 0.05 /iMpc^^: this 
is typically the scale where it is assumed to be valid. This casts 
doubts about the reliability of the works in which linear the- 
ory predictions are used to constrain modified gravity theories 
such as chameleon, dilaton, symmetron and f{R) gravity. 

Another feature of the chameleon theory is its efficiency of 
screening. The model parameters here, such as /3o and ^, are 
chosen to be roughly the same as those in our previous dila- 
ton and symmetron simulations |fl2|], but whilst the nonlinear 
matter power spectra in those models can differ from those in 
ACDM by more than 30 — 40%, chameleon theories generally 
predict much smaller deviations (< 10%), indicating that the 
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FIG. 10. The same as Fig|9l but at z = 1. 



chameleon screening could restore GR much more easily. For 
the same reason, the effect of the fifth force also diminishes 
more quickly backwards in time, compared to the symmetron 
and dilaton cases lfl2ll - indeed at redshift z ~ 1 the fifth force 
is almost completely screened in all our simulated models ex- 
cept for Bl, which has r = 2.0, meaning that the scalar field 
mass m increases more slowly with matter density. The result 
implies that the strength of the fifth force is very sensitive to 
r, which is, of course, as expected. 

Similar features can also be seen from the dark matter halo 
mass functions. Here we find that, compared with the dila- 
ton and symmetron theories lfl2ll . the deviations from ACDM 
are more suppressed in the high-mass end, which can be be- 
cause large halos are more efficient in self-screening and also 
tend to be more screened by the environment because they are 
more likely to live in high-density environments. This is qual- 
itatively similar to what we see in f{R) gravity simulations 
ifsill . Notice that the time evolution of the halo mass function 
shows the same pattern as the nonlinear matter power spectra, 
namely that at z = 1 the deviation from ACDM is very small. 

The high efficiency in chameleon screening means that our 
choices of the parameter values might be too conservative: a 
deviation from the ACDM matter power spectrum of < 10% 



can hardly be detected with precision in the near future, espe- 
cially because the deviations are mostly on small scales where 
baryonic physics and other effects could already be important. 
Consequently, we think that future simulations of chameleon- 
type theories should be done for less conservative choices of 
parameters, namely larger values of /3o , s, ^ and smaller values 
for r. We hope that this work can serve as a useful guidance 
for such future works. 
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